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The flow in a cylinder driven by time harmonic oscillations of the rotation rate, 
called longitudinal librations, is investigated. Using a theoretical approach and ax- 
isymmetric numerical simulations, we study two distinct phenomena appearing in this 
librating flow. First, we investigate the occurrence of a centrifugal instability near 
the oscillating boundary, leading to the so-called Taylor- Gortler vortices. A viscous 
stability criterion is derived and compared to numerical results obtained for various 
libration frequencies and Ekman numbers. The strongly nonlinear regime well above 
the instability threshold is also documented. We show that a new mechanism of 
spontaneous generation of inertial waves in the bulk could exist when the sidewall 
boundary layer becomes turbulent. Then, we analyse the librating flow below the 
instability threshold and characterize the mean zonal flow correction induced by the 
nonlinear interaction of the boundary layer flow with itself. In the frequency regime 
where inertial modes are not excited, we show that the mean flow correction in the 
bulk is a uniform rotation, independent of the Ekman number and cylinder aspect 
ratio, in perfect agreement with the analytical results of Wang [J. Fluid. Mech., 41, 
pp. 581 - 592, 1970]. When inertial modes are resonantly excited, the mean flow 
correction is found to have a more complex structure. Its amplitude still scales as 
the square of the libration amplitude but now depends on the Ekman number. 
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I. INTRODUCTION 



Rotating flows generically support oscillatory motions called inertial waves 1 ' 2 , whose fre- 
quencies range between plus and minus twice the spin frequency. When the rotating flow 
is limited by boundaries, these waves can combine to form inertial modes whose structure 
and properties depend on the geometry of the container. These inertial modes are well- 
characterized in simple geometries only, as the cylinder or the sphere 1 . In other geometries, 
as the spherical shell for instance, the singular structure that some of them exhibit is still 
the subject of active research 3-5 . Usually damped by viscosity, inertial modes can never- 
theless be excited by various harmonic forcings such as the natural processes of libration, 
precession and tides in planetary fluid cores, providing that their azimuthal period m and 
temporal frequency u are close to those of the forcing. Once a mode is resonantly excited, 
its nonlinear self-interaction generates an intense axisymmetric flow referred to as the zonal 
flow (see Morize et al. 6 ). In the absence of inertial mode resonance, harmonic forcing can 
also generate a zonal flow through nonlinear effects in the Ekman layer 7-14 . The study of the 
zonal flow with and without resonance is of interest in a geophysical context, as it provides 
constraints on the possible dynamics of liquid cores, internal oceans and atmospheres of 
planets 15 . It is also important in the aeronautical context of rotating flying objects 16 since 
the internal fluid dynamics could influence their trajectory. 

In this paper, we focus on a particular forcing: the time harmonic oscillations of the 
rotation rate of a container, called longitudinal librations. The early studies of librational 
forcing mainly focused on the excitation of inertial modes in a sphere or spherical shell. 
The first experimental study of Aldridge & Toomre 17 confirmed the theoretical resonance 
of inertial modes in a sphere and evidenced the presence of Taylor-Gortler vortices near the 
outer boundary with axes oriented in the azimuthal direction of the sphere 18 . These exper- 
imental results were then confirmed numerically by Rieutord 19 and extended by Tilgner 20 
who studied the linear response to a librational forcing in a spherical shell and investigated 
the attractors associated with these inertial modes. More recently, Noir et al. 21 investigated 
experimentally by direct flow visualization, using Kalliroscope particles, the appearance of 
Taylor-Gortler vortices in a spherical container. This study was then complemented by the 
numerical work of Calkins et al. 9 who also found a nontrivial zonal flow in the interior. 
By the same time, using methods developed for a precessing sphere 12 , Busse 7 obtained an 
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analytical solution of the flow in a librating sphere in the limit of small libration frequency, 
when inertial modes are not excited. This weakly non-linear theory was then confirmed by 
Sauret et al. 8 using a combined numerical and experimental study. 

In addition to these studies in spherical geometries of direct geo/astrophysical interest, 
the case of a librating cylinder has also been investigated. In a little known paper, Wang 11 
already demonstrated theoretically and experimentally that the time-harmonic forcing gen- 
erates a zonal flow in the bulk. He showed that this flow comes from nonlinear interactions 
in the top and bottom Ekman layers of the cylinder. He also provided an estimate for the 
amplitude of the zonal flow in the interior for arbitrary values of the frequency of libration in 
the absence of inertial mode forcing. The limit of small libration frequency compared to the 
mean rotation rate was recently recalculated by Busse 14 . Noir et al. 10 studied experimen- 
tally, via direct flow visualization and LDV measurements, the mean zonal flow at libration 
frequencies comparable to spin rate. They also analysed the generation of Taylor-Gortler 
vortices and suggested a scaling law for the appearance of the vortices different from that 
in the sphere. More recently, Lopez & Marques 22 performed a three dimensional numerical 
study of a librating cylinder at moderate Ekman number (E = 1CT 4 ) for three values of 
libration frequency and various values of the libration amplitude. They observed the spon- 
taneous generation of inertial modes for libration frequency above u > 2 and proposed a 
mechanism for the excitation of inertial modes based on a period doubling process. The 
present paper comes as a complement to these works. 

The paper is organized as follows. In section 2, we introduce the governing equations, the 
notations and the numerical method. In section 3, we study the occurrence of Taylor-Gortler 
vortices at the outer boundary. We derive a local stability criterion which is systematically 
tested as a function of the libration amplitude, the libration frequency and the Ekman 
number. Then we compare this criterion to numerical results and available experimental 
data. The spontaneous generation of inertial waves induced by the Taylor-Gortler vortices 
is discussed in section 4. Section 5 is devoted to the analysis of the zonal flow over a wide 
range of libration frequencies, u G [0.05; 25]. We show that in the absence of inertial mode 
resonance, the numerical results are well described by the analytical prediction of Wang 11 
in the bulk. Configurations where inertial modes are present are also documented in this 
section. 
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II. MATHEMATICAL FORMULATION AND NUMERICAL MODEL 



A. Governing equations 

We consider a homogeneous and incompressible fluid of kinematic viscosity v in a cylinder 
of radius R and height H (see figure 1) rotating at the mean rotation rate f2 along its 
symmetry axis. Using R and VLq^ 1 as the length scale and the time scale, respectively, the 
rotation rate of the cylinder is given by 

Q{t) = [l + e cos(ut)]e z , (1) 

where u = uub/Qo is the dimensionless libration frequency, e = A0u is the amplitude of 
the librational forcing with A6 the libration angular amplitude, and e z is the unit vector in 
the direction of the rotation axis. The resulting dimensionless period of libration is given 
by T Ub = 2ir/u. In the frame rotating at the mean rotation rate, the equations governing 
the fluid motion are 

u 

— + (u- V)u + 2e z x u = -VU + EV 2 u, (2) 

at 

V • u = 0, (3) 

where E is the Ekman number defined by 

(4) 

u is the velocity of the fluid in the rotating frame and II is the reduced pressure taking into 
account the centrifugal force. The fluid satisfies no-slip boundary conditions on the cylinder 
walls: 

u = e cos(ut) e e at r = 1, (5) 
u = er cos(wt) e e at z = ±a/2, (6) 

where a = H/R is the aspect ratio of the cylinder, (r, z) being the cylindrical radial and 
vertical coordinates. 



In this study we consider the limit of small Ekman number, E 1, and we assume that 
spin- up does not occur in the interior during a libration cycle, i.e. u ^> \fE: this means 
that in the bulk, at first order, the fluid remains in solid body rotation and does not follow 
the librating outer boundary. 
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FIG. 1. Schematic view of the cylindrical cavity and the cylindrico-polar coordinate (r,6,z). The 
z-axis is the symmetry axis and of rotation. The meshgrid used for the numerical simulations is 
also represented, with a denser mesh near the axis of rotation and near the outer boundaries (see 
section II B for details) . 

B. Numerical model 

Equations (2) and (3) with the boundary conditions (5) - (6) are solved with the com- 
mercial software, Comsol Multiphysics®. Axisymmetric simulations are performed using a 
standard finite element method in which the numerical grid is composed of two subdomains 
(see figure 1): (i) a boundary layer domain typically composed of 8 layers of quadriangular 
elements all along the outer boundaries and around the axisymmetric axis; (ii) a bulk zone 
with triangular elements. 

All elements are of standard Lagrange PI — P2 type (i.e. linear for the pressure field 
and quadratic for the velocity field). Note that the finite element method does not induce 
any particular problem around r = and that no stabilization technique has been used in 
this work. The temporal solver is the so-called Implicit Differential-Algebraic solver (IDA 
solver), based on backward differenciating formulas 23 . 

At each time step the system is solved with the sparse direct linear solver PARDISO 24 . 

Starting at time t — with a fluid at rest in the rotating frame, libration of the outer 
boundary is turned on and computations are pursued during a few spin- up times, E^ 1 ^ 2 , 
until a stationary state is obtained. We then determine the zonal flow ztg by averaging over 
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FIG. 2. Convergence with the number of Degrees of Freedom (DoF) of the averaged rotation rate 
of the zonal flow in the interior, rescaled by e 2 (see section V). The simulations are performed for 
E = 5.1CT 5 , e = 0.1, u = 0.1 and an aspect ratio a = 2. The black dashed line corresponds to 
113 387 DoF which is used in the following. The value < uq > is obtained by averaging over three 
periods the azimuthal velocity field and then taking the mean value of the zonal flow between r = 
and r = 0.5 at z = 0.2 which leads to —0.1496 for the highest resolution. 

three libration periods the numerical velocity field at a fixed axial location. The numerical 
model has been previously benchmarked and used in Sauret et al. 8 to study the steady 
zonal flow induced in a librating sphere in the limit of low libration frequency. The number 
of degrees of freedom (DoF) used in the simulations is 113 387 DoF which is sufficient to 
ensure accurate results (see the convergence test in figure 2) while keeping a reasonable CPU 
time (i.e. about 8 hours on a standard workstation). Using this method, we systematically 
explored the range 0.05 < u < 25 with 0.01 < e < 0.85 and E as small as E = 2 ■ 10~ 5 . 




III. CENTRIFUGAL INSTABILITY 

As already documented 10 for fixed values of u and E, different flow regimes are expected 
near the sidewall boundary depending on the forcing amplitude e. If e is smaller than a 
critical amplitude e c , the flow remains stable in the whole cylinder (figures 3(a) and 3(b)). 
Then, increasing the amplitude e beyond the critical value e c leads to the growth of an 
instability near the outer boundary (figures 3(c) and 3(d)). This instability corresponds to a 
centrifugal destabilization of the time-periodic azimuthal velocity field in the boundary- layer 
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FIG. 3. Visualizations of the libration flow in the upper half of the cylinder (0 < r < 1; < z < 1): 
axial velocity u z (upper pictures) and azimuthal velocity uq (lower pictures) at the same time 
during the retrograde phase (t = 0.4015 All simulations have been done for a = 2, E = 5.10 -5 
and u> = 0.05. (a) and (b) show a stable case where no instability is observed, (c) and (d) are 
centrifugally unstable with the generation of Taylor-Gortler vortices, (e) and (f ) illustrate the more 
complex flow that takes place for a larger excitation amplitude e. 

generated by libration. It is characterized by the appearance of counter-rotating Taylor- 
Gortler vortices close to the sidewall boundaries. For larger values of e, the flow becomes 
more complex (figures 3(e) and 3(f)). We have not tried to fully characterize this last flow 
regime because non-axisymmetric features inaccessible to our axisymmetric simulation are 
then expected to be present (but see section IV). 

Figure 4(a) and (b) show typical temporal signals of the axial velocity near the sidewall 
boundary in stable and unstable regimes. In both figures, the azimuthal velocity of the 
sidewall is also plotted. In figure 4(a), the axial velocity is taken at the location r = 0.99, 
z = that is very close to the boundary and on the horizontal plane of symmetry of the 
cylinder. Below the instability threshold (for e = 0.28), no signal is observed at this point 
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FIG. 4. Axial velocity u z for e = 0.28 (blue dashed-dotted line, left scale), e = 0.34 (red continuous 
line, left scale) and azimuthal velocity u s id e of the sidewall in the rotating frame (black dashed 
line, right scale) in a cylinder of aspect ratio a = 2 (a) at the point (r = 0.99; z = 0) (b) at the 
point (r = 0.99; z = 0.1). Here, E = 4.10" 5 and u = 1. 

whereas a signal oscillation at the libration period is observed at z — 0.1 (figure 4(b)). These 
oscillations, which are visible when we are not in the plane of symmetry, come from the 
Ekman pumping from top and bottom boundary layers during the retrograde and prograde 
phases. 

When e is increased above the instability threshold, the axial velocity signal is modified. 
Either it becomes visible (at z = 0; figure 4(a)) or it reaches larger amplitude (at z = 0.1; 
figure 4(b)). It can be associated with the growth of the centrifugal instability at the end 
of the retrograde phase. A net peak is indeed observed in the middle of the prograde phase 
at both locations. 

In a previous study, using direct flow visualization with Kalliroscope particles, Noir et 
al. 21 defined a boundary layer Reynolds number Rcbl & s the relative strength of inertial 
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and viscous forces in the boundary layer: ReBL = u5/u where v is the kinematic viscosity 
of the fluid, u is the characteristic azimuthal velocity and 5 is the boundary layer thickness. 
For Re BL > RtBL,cr where ResL,cr is a critical boundary layer Reynolds number obtained 
experimentally, a centrifugal instability develops near the sidewall. Using scaling arguments, 
they established that ReBL,cr should scale as eE~ x l 2 in the spherical case, which was then 
confirmed numerically by Calkins et al. 9 . The influence of the libration frequency uj was 
not studied. In the cylinder, based on laboratory experiments, Noir et al. 10 suggested that 
the critical boundary layer Reynolds number should scale as e£~ 3 / 4 . In the following, a 
theoretical expression for the onset of Taylor-Gortler vortices is derived using local stability 
considerations. This expression takes into account the dependence with the frequency and 
reconciliates the cylindrical and spherical geometries. 



A. Theoretical determination of the threshold of the centrifugal instability 

Let us consider the situation presented in section II. Assuming that the amplitude of 
excitation is small, e 1, in the rotating frame, one can expand the velocity field in power 
of e. At the order e, equation (2) writes 

C ^- + 2e z xu 1 = -Vp 1 + E W, (7) 

where the subscript 1 denotes the terms of first-order in e. By performing a classical bound- 
ary layer analysis 11 , we obtain the azimuthal correction of the velocity field in the sidewall 
boundary layer: 

ug 1 = eexp ( -=- t — J cos I out -=L- J — ] . (8) 




VE V 2 

This boundary layer flow of thickness 0(y/E/u) corresponds to what is called a Stokes 
layer. As long as the period of oscillation is small compared to the spin-up time (that is 
u ^> V~E), this boundary layer flow does not generate any flow in the bulk of the fluid 11 ' 14 . 
By contrast, the boundary layer flow on top and bottom boundaries induces a E 1 ^ 2 correction 
of the flow in both the bulk and the sidewall boundary layer (see Wang 11 ). This weak axial 
pumping flow remains small when cu ^> \[E and is neglected in the present analysis. 

The experimental results 10 tend to demonstrate that the flow (8) can become centrifugally 
unstable with respect to axisymmetric perturbations. Is is then natural to apply the Rayleigh 



criterion, which is a necessary 25 and sufficient 26 condition for instability of axisymmetric 
inviscid perturbations in a stationary flow. This criterion is 



^ d / 2 (tot) 2 \ 



<0, 



(9) 



where u e ° is the total azimuthal velocity field in the absolute frame of reference. This 
criterion can be satisfied only in the boundary layer where the Rayleigh discriminant $ 
reduces at first order to 



By using (8) for u\ in (10), we implicitly perform a quasi-steady analysis. In other words, 
we assume that the time variations of the flow are slow compared to the temporal scale 
associated with the centrifugal instability. The instability criterion is thus local in time. 
By considering only the moment where the base flow is the most unstable, we then expect 
to provide a lower boundary for the instability domain. Here, the most unstable profile is 
reached at t m — 3 n/ (4 a;) and the minimum value of the Rayleigh discriminant $ is obtained 
at r rn = 1 and equals: 



When & m in > 0, the Rayleigh criterion is not satisfied and we thus do not expect insta- 
bilities. In terms of e, this means that if 



the flow is stable. This criterion is non-viscous and does not take into account the viscous 
damping that will affect the small scale instability modes. To obtain a viscous criterion, it is 
necessary to solve the viscous stability equation for the flow given by (8). Taking into account 
that the most unstable profile is obtained at t m , a lower limit for the stability threshold can 
be obtained by analysing the stability properties of such a profile. Assuming both e <C 1 
and y/E/u < 1 with e^/u/E = 0(1), the viscous stability problem for axisymmetric 




(10) 




(11) 




(12) 
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perturbations (u r , u$ } u Zj p) exp (i k z + a t) reduces to solving the system: 

<TUr-2u e =-^ + - [-Q^-k 2 \ u r , (13) 



™o+ \u r = - ^--k 2 u e , (14) 



_ dul(r)\ uj ( d 2 

~dT) Ur= 2 \dr 2 

uj ( d 2 \ 
au z = -ikp+- I — -k 2 \ u g , (15) 

du r , . . 

—^ + 1 ku z = 0, (16) 

with the boundary conditions 

u r = ug = u z = at f = 0,+oo, (17) 



where 



1 — r [uj 



E V 2' 

z [uj 



\[E V 2 



(18) 
(19) 



P=^Jl (20) 



£ V 2' 

; (UJ 

IV 2' 



*=-W£ (2D 



uj(f) = e _fi v / 2cos(f). (22) 

This eigenvalue problem for a is solved numerically using a pseudospectral method. The 
radial discretization is based on Chebyshev functions. The parameter E has disappeared 
from the equations. This means that the eigenvalues a obtained by solving (13)-(16) with 
(17) are function of k, e and uj only. This implies that the maximum growth rate over all 
possible A; is a function of e and uj only: 

a max = max a = a max (e, u) (23) 

k 

The stability curve, defined by a max = 0, is then given by a curve e = e c (oj) which is a 
function of uo only. Using the primitive variable, the viscous criterion for stability therefore 
reads 

e < & = ~ €c{u}) (24) 
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FIG. 5. Critical amplitude e c /y/E for the centrifugal instability obtained by the theory. Solid line: 
viscous criterion e^\uj)/y/E; dashed line: inviscid criterion ei NV \uj)/\^E = 2/y/oj 



In figure 5, we have plotted e^{u)/yfE as a function of uj. This plot is compared to the 
non viscous prediction ei NV \u>)/y/E = 2/y/uj. Note that when uj — > 0, e£ ~ ei N . This 
is in agreement with the system of equation (13)-(16) in which the viscous terms vanish as 
u — >■ 0. However, both ei NV ^ and 6c diverge as u — > 0. This divergence is associated with 
the blow up of boundary layer thickness when oj — > 0. In fact, this boundary layer does not 
blow up but transforms into a Stewartson layer of thickness E l l 4 when to becomes of order 
E 1 / 2 or smaller. Under the hypothesis u ^> \[E we never reach this regime. Moreover, 
we see that neglecting the viscosity leads to a non-viscous criterion which decreases while 
increasing w, in contradiction with the experiments (see below). On the contrary, the viscous 
criterion has a correct tendency to increase with u. 

The above viscous stability criterion has been obtained by considering the most unstable 
profile in a frozen state. It should then, a priori, provide a lower limit for the critical desta- 
bilization amplitude. In principle, as the base flow oscillates, a Floquet stability analysis 
should be performed to obtain a more precise estimate for the critical amplitude but we 
discuss below why we do not think that such a development is worthwhile. Finally, note 
that the above viscous stability criterion has only assumed that the boundary layer is lo- 
cally vertical. It should therefore apply close to the vertical boundary of any axisymmetric 
containers. In particular, we expect the same criterion to apply for the flow in a sphere 
near the equator. In the following sections, the viscous stability criterion is compared to the 
numerical simulations and experimental observations. 
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B. Influence of the Ekman number 

Our numerical simulations allow us to recover the range of parameters explored in previ- 
ous experiments in the literature 10 but also to increase the range of studied Ekman numbers. 
Moreover we can explore more precisely the transition between the stable/unstable regimes 
with the measurement of the axial velocity. Finally, we can study the influence of the li- 
bration frequency which was not investigated previously. The method used to determine 
the stability of the flow is the same as the one used by Calkins et al. 9 : for given E and u> 
we start at an arbitrary amplitude of libration e and observe if the flow is stable (respec- 
tively unstable), then we increase (respectively decrease) the amplitude of libration until we 
reach an unstable (respectively stable) regime. The determination of the instability is done 
both by direct visualization, and by plotting the vertical velocity near the outer boundary 
at some locations as it has been observed that an important variation of its amplitude is 
visible at the transition between stable and centrifugally unstable regimes. If the flow is 
stable, the vertical velocity remains small and its amplitude is periodic with a period equal 
to the libration period (see figure 4(a)). Once we reach the unstable regime the vertical 
velocity suddenly increases (figure 6) and becomes less regular (see figure 4(b)). The critical 
libration amplitude e c is defined as the largest stable e. We have not observed any hysteresis. 
We therefore expect the instability to be supercritical as for classical Taylor-Couette flow. 

Figure 6 shows the maximum peak to peak amplitude, i.e. the maximum value minus 
the minimum value, of the vertical velocity as a function of the amplitude of excitation e 
for a given Ekman number and a given libration frequency. Due to the presence of Ekman 
pumping from top and bottom boundaries, even below threshold, the value of the axial 
velocity depends on the location of the measurement. The axial velocity is weak close to the 
plane of symmetry but large close to the top and bottom boundaries. A transition is visible 
close to e c = 0.355. This transition is very clear in the measurements at z = 0.1 as the 
axial velocity increases by a factor three. Visualizations also show that the unstable rolls 
do not always develop at the same location in the cylinder. For small libration frequencies, 
typically u> < 1, we observe that the Taylor- Gortler vortices first develop near the equator of 
the cylinder, i.e. at z ~ 0. Increasing the frequency of libration leads to the first development 
of these unstable rolls at mid-distance between the corner and the middle plane. This feature 
is probably also related to the Ekman pumping from top and bottom boundaries. 
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FIG. 6. Amplitude of the axial peak-to-peak velocity (Au z ) max in a cylinder of aspect ratio a = 2 
for E = 6.10" 5 , oj = 0.8 and various values of e. The measurements are performed at different 
locations: the blue circles at (r = 0.95; z = 0.1), the red squares at (r = 0.95; z = 0.4) and the 
green triangles at (r = 0.95; z = 0.9). We see a clear transition between stable/unstable regimes 
at e c = 0.355. 

Figure 7 presents the results of our numerical investigation of the Ekman number depen- 
dence. The £' 1//2 -scaling law is clearly observed and is well fitted by the viscous stability 
threshold (24) which gives e c = 47.8 E 1 ! 2 at u — 1 with no adjusting parameter. If we take 
5 = \fE^j2Vt a /ujR as spatial scale and u = y/2eil R/2 as velocity scale, as prescribed by 
the boundary layer theory (see (8)), this corresponds to a critical boundary layer Reynolds 
number ReBL,cr = 47.8. This is surprisingly close to the value RtBL,cr = 46 obtained in a 
spherical shell by Calkins et al. 9 at the same u. 

In figure 8, we have reported the experimental observations of Noir et al. 21 for the same 
aspect ratio as in figure 7. The theoretical prediction and the scaling law suggested by 
Noir et al. 10 are also shown. The agreement with the theory is not as good as with our 
numerical simulations. In particular we do not understand the fact that the experimental 
threshold is well below our estimate which is expected to be a lower limit. The same 
observation was made in the sphere between experimental 21 and numerical 9 results. The 
discrepancy between these measurements was explained by Calkins et al. 9 by a systematic 
higher noise level present in the experimental setup. It could also be due to the visualization 
with Kalliroscope particles which do not directly measure the velocity field. Note that a 
scaling in E 1 ! 2 still fits the experimental data with a good agreement. Therefore, regarding 
the experimental measurement both scalings in E 1 / 2 or E 3 ^ could a priori be possible. 
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FIG. 7. Stability diagram for the onset of Taylor-Gortler vortices for a libration frequency of oj = 1 
in a cylinder of aspect ratio a = 2. The black bars are the numerical results together with the 
uncertainty in the estimation of the threshold. The blue continuous line is the viscous stability 
criterion (24) without adjusting parameter: 47.8 E 1 / 2 and the red dashed line is the best fit line 
using the scaling of Noir et al. 10 : 543.2 E 3 / 4 . 

However, the latter scaling in E 3 / 4 does not hold for the numerical results (see figure 7). 



C. Influence of the libration frequency u 

The influence of the libration frequency u on the instability threshold has not been 
considered in previous studies. Here we consider 8 different values of u arbitrarily chosen 
in the range u G [0;2]: 0.1, 0.2, 0.3, 0.4, 0.65, 0.8, 1 and 1.5. The results are reported in 
figure 9. A clear increase of the instability threshold with u is observed in the simulation. In 
this figure, we have also plotted the viscous stability criterion (24). Whereas the theoretical 
prediction captures correctly the global increase of the threshold with u, it fails in providing 
correct quantitative values for small and large u. 

For small u, we predict a threshold which is twice smaller than the numerically observed 
value. We cannot attribute this difference to the time variation of the base flow, which was 
not taken into account in the theory. Indeed for small u, the quasi- st at ionnary approximation 
applies and we expect a Floquet stability analysis to provide a similar stability threshold. 
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FIG. 8. Experimental results from Noir et al. 10 obtained by direct visualization plotted in terms 
of the libration forcing e versus E for a fixed value of ui = 1. The blue squares are centrifugally 
stable and the red circles are centrifugally unstable. The black line is the viscous stability criterion 
(24), the blue dashed-dotted line is the best fit line in E 1 / 2 : 22 E 1 / 2 , and the green dashed line is 
the fit line suggested by Noir et al. 10 : 250 E 3 ^. 
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FIG. 9. Stability diagram for the onset of Taylor- Gortler vortices as a function of the frequency of 
libration u. The squares are the numerical simulations for E = 6.10 -5 in a cylinder of aspect ratio 
a = 2 and the dashed line is the viscous stability criterion (24) with no adjustment parameter. 



Instead, we suspect the difference to be associated with a bad description of the boundary 
layer flow. For small frequencies (if the condition u ^> \/E is no longer satisfied), transient 
Stewartson layers generated by top and bottom boundaries could have time to appear 27 and 
modify the Stokes layer solution (8). 
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For large u, our theoretical prediction is above the numerical stability threshold. Again, 
a more precise Floquet analysis would have therefore not explained the difference. Indeed, 
the global stability results obtained by a Floquet analysis would have provided a smaller 
instability domain and therefore an even larger value for the threshold. Here we think that 
the discrepencies come from the nonlinear effects. Our theoretical approach which assumes 
a linear boundary layer flow is indeed expected to break down when e is not sufficiently 
small. 

IV. SPONTANEOUS GENERATION OF INERTIAL MODES 

In the previous section, we have studied the occurrence of the centrifugal instability on the 
sidewalls and have mainly considered the instability threshold. In this section, we consider 
configurations well above the threshold for which rich dynamics are observed. Our numerical 
simulations tend to demonstrate that increasing the libration amplitude e eventually leads to 
a turbulent flow in the boundary layer. We have also observed that these turbulent regimes 
can drive motion in the bulk under the form of spontaneously generated inertial waves. 

In figure 10, we have plotted d z u r for configurations above the instability threshold for 
different frequencies. In all cases, we are well above the threshold, i.e. in regimes where the 
boundary layer can be considered as turbulent. We observe structures in the bulk of the 
fluid, which seem to share some common features. In particular, we observe in all cases, 
small scale structures oriented at a fixed angle ~ 45 ° . We think that these structures are 
inertial waves excited by the turbulent motion in the boundary layer. In figure 11, we have 
plotted a typical frequency spectrum obtained for each case. For each configuration, we do 
observe a peak at the frequency of excitation, that is u = 2.5 in (a), u = 0.1 in (b), and 
oo = 4.1 in (c). Each case also possesses a peak at approximatively cu « 1.4 which corresponds 
to the frequency of the inertial wave pattern oriented at 45°. Indeed, the orientation pattern 
is also the direction of the group velocity of the inertial wave packet which is related to its 
frequency by the relation 1 u = 2sin#. This relation gives u ~ 1.41 for 9 ps 45°. 

The case (a) of a forcing frequency ui = 2.5 was previously studied by Lopez & Marques 22 . 
They also found coherent structures in the bulk. They argued that those structures are 
inertial modes oscillating at ui/2 which are excited by a period-doubling mechanism. Here 
this mechanism gives a frequency u ex = 1.25 which corresponds to the second peak in figure 
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(a) (b) (c) 

FIG. 10. Axial derivative of the radial velocity in the librating cylinder of aspect ratio a = 2 at 
time t = 3vr/(2o;) for (a) uj = 2.5, e = 0.70, (b) u = 0.1, e = 0.6, (c) uj = 4.1, e = 0.85. All 
simulations are performed for E = 3.10 -5 . The amplitude of the colorbar is chosen to visualize 
inertial modes (a saturation of the axial velocity near the outer boundary of the container is thus 
visible) . Dotted lines correspond to an orientation of 45 ° . 



11(a). The pattern observed in figure 10(a), especially the two inertial wave beams emitted 
from the corners, does not correspond to this inertial mode. We think that it is associated 
with a stronger turbulence activity near the corners for this particular forcing frequency. 

The case (b) of a forcing frequency uj — 0.1 corresponds to the situation where the 
inertial wave pattern is the most visible and the most regular. The frequency spectrum also 
possesses a larger peak at uj w 1.4 in that case. Note that other peaks (at uj = 0.2,0.3, ..) 
corresponding to harmonic frequencies of the forcing are also visible on the spectrum. 

The case (c) of a forcing frequency uj = 4.1 is interesting as the period doubling mechanism 
proposed by Lopez & Marques can not be active since uj/ 2 is not in the inertial frequency 
range (0,2). The inertial wave pattern can not therefore be directly related to the forcing 
frequency. As for the two other frequencies, a pattern preferentially oriented at 45° is still 
observed. 
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FIG. 11. Frequency spectrum in a cylinder of aspect ratio a = 2 for E = 3.10 -5 . The three 
cases corresponds to the situation presented in figure 10, namely for (a) to = 2.5, (b) to = 0.1, (c) 
to = 4.1. The spectrum is obtained by axially averaging the radial velocity at r = 0.6. The black 
dotted line indicates in each cases the libration frequency 



The three different cases have been simulated for different libration frequencies, to = 0.1, 
to = 2.5 and to = 4.1. Independently of the frequency of the forcing, the inertial waves 
generated by the turbulent flow near the outer boundary always exhibit a peak in the 
frequency spectrum around to ~ 1.4. It suggests that another mechanism than direct forcing 
or period doubling could exist for the generation of inertial waves in a cylinder. Because the 
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characteristic of the inertial waves do not depend on the forcing frequency, we suspect that 
the generation could be due to the turbulent motion in the boundary layer. Such a forcing 
by turbulence has not been studied in detail in the context of rotating flows. Nevertheless, 
we can reasonably think that the results obtained in the context of gravity waves where 
similar observations were made could apply (see Dohan & Sutherland 28 ). More works on 
this interesting issue which is beyond the scope of the present paper are clearly needed. 



V. MEAN ZONAL FLOW INDUCED BY THE LIBRATIONAL FORCING 



In this section we analyse the mean zonal flow induced by the librational forcing that is the 
time averaged value of the azimuthal velocity. We first consider configurations where inertial 
modes are not excited in order to compare our numerical results to a previous analytical 
solution given by Wang 11 . Then we analyse the mean zonal flow when inertial modes are 
present in the interior and describe the key features of the resulting flow in this regime. We 
assume that the amplitude of libration is sufficiently small such that the sidewall boundary 
layer remains centrifugally stable. 

A. Zonal flow for w < 1 or w > 2 

A full analytical solution of the zonal flow in a librating cylinder has been obtained by 
Wang 11 for an arbitrary libration frequency u, assuming that no inertial mode is excited 
in the interior. He showed that the mean flow correction is at leading order an azimuthal 
flow of the form ug = e 2 rfi 2 (w), where Q 2 is uniform in the bulk of the fluid. The explicit 
expression of Q 2 in terms of the frequency is given in the appendix. 

Wang 11 predicted that this mean flow should match the evolution of the sidewall boundary 
where Ug = in a E 1 ^ viscous layer. 

Note that Wang 11 also demonstrated that a E 1 ^ 3 viscous layer is needed to satisfy the 
continuity of the axial velocity field. A composite approximation for the azimuthal velocity, 
which takes into account the E l / A layer solution, is obtained in the form 



It is this expression, with Vt 2 given by (27) and (28) from the appendix, that we shall test 
below. 




(25) 
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FIG. 12. Angular velocity 0,2 = ue/( r ^ 2 ) of the zonal flow in the whole cylinder of aspect ratio 
a = 2 obtained numerically for e = 0.05, E = 4.10 -5 , (a) uj = 0.05 and (b) u = 6. The zonal flow 
is obtained by averaging the azimuthal velocity over 3 libration periods. 

In figure 12, the mean angular velocity obtained from the numerical simulation is plotted 
in a meridional plane (r, z) for two libration frequencies. We clearly see in these plots that 
the rotation is uniform in the bulk in both cases. Moreover, near the sidewall (r = 1), 
the mean angular velocity is primarilly axially independent on z, in qualitative agreement 
with formula (25). Note however that the small localized structures in the corners are not 
described by Wang's theory. 



Figure 13 shows the dependence of the mean angular velocity with respect to e: all 
the numerical curves for VI2 = ue/ (r e 2 ) superimpose as e is varied which demonstrates the 
scaling in e 2 of the mean zonal flow. In figure 14, we have plotted the rescaled angular 
velocity Q2 = ug/(r e 2 ) as a function of p = (1 — r)/E l / A keeping e and to fixed but varying 
the Ekman number for two libration frequencies. We see that for large p (that is far from 
the sidewall) all the curves collapse approximatively to the theoretical prediction. This 
means that the mean angular velocity in the bulk does not depend on the Ekman number 
as predicted by the theory. Close to the wall, discrepancies between the theory and the 
numerics are by contrast visible: in the absence of instabilities, a prograde peak whose 
amplitude varies with the Ekman number is obtained. We suspect that the small scale 
structures observed near the corners are at the origin of these discrepancies. Indeed, Wang 11 
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FIG. 13. Angular velocity 0,2 = ue/( r ^ 2 ) of the zonal flow as a function of the radius r for a 
cylinder of aspect ratio a = 2 and E = 5.1CP 5 , uj = 0.05 and various amplitudes of libration e: 
e = 0.01 (cyan dotted line), e = 0.05 (red dashed line), e = 0.1 (green dashed-dotted line), e = 0.2 
(blue continuous line). The black thick continuous line shows the theoretical expression of Wang 11 
including the £' 1 / 4 -layer solution 11 [expression (25)]. 

did not solve the corner region where sidewall and top/bottom boundary layers encounter. 
He assumed that both layers could be treated independently. The present numerical results 
tend to show that this hypothesis could be incorrect and that the sidewall region could be 
strongly affected by the corners. Characterizing the flow close to the corner remains a great 
challenge which is beyond the scope of the present paper. 

The value of the mean angular velocity obtained in the bulk has been systematically com- 
puted in different locations for frequencies ranging from 0.05 to 25. The results are reported 
and compared to the theory in figure 15. When u > 2, very good agreement is obtained 
between theory and numerics. In this regime, no inertial modes can be directly excited by 
the forcing. The hypothesis of Wang (i.e. no inertial mode forcing) is then satisfied. Note 
in particular that the zonal flow in this limit can be prograde. Good agreement is also 
observed for the smallest frequency (u = 0.05) that we have considered. For this frequency, 
we suspect that no inertial mode is excited due to the finite length of the cylinder and the 
relatively large value of our Ekman number. By contrast, for the other frequencies between 
and 2, inertial modes are present in the system which lead to differences between theory 
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FIG. 14. Angular velocity Q2 = u$/(re 2 ) of the zonal flow as a function of the local coordinate 
p = (1 — r)/.E 1//4 in a cylinder of aspect ratio a = 2 plotted for e = 0.1, (a) w = 0.07 and 
(b) lj = 3 and various Ekman number E: E = 10~ 4 (green dashed-dotted line), E = 8.10 -5 
(blue continuous line), E = 4.10 -5 (red dashed line). The black thick continuous line shows the 
theoretical expression of Wang 11 including the i? 1//4 -layer solution 11 [expression (25)]. 

and numerical results (see figure 15(b)). This regime in which inertial modes are excited is 
studied in the next section. 



B. Zonal flow in the presence of inertial modes 

In this section, we shall try to provide some generic features of the zonal flow in the 
presence of inertial modes. Again, we consider small amplitude configurations where the 
boundary layer remains centrifugally stable. 
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FIG. 15. (a) Angular velocity = W/( re2 ) of the zonal flow, as a function of the frequency of 
libration in a cylinder of aspect ratio a = 2 for E = 4.10 -5 and e = 0.05. The black continuous 
line is Wang's prediction [expression (27) or (28)], the blue squares are the numerical results at 
(r = 0.2; z = 0.2), the red circles at (r = 0.4; z = 0.2) and the green stars at (r = 0.6; z = 0.2). 
The discrepancy for uj £ [0.1; 2] is likely due to the contribution of inertial modes, (b) is a zoom of 
(a) within the region were inertial modes are expected. 

Experimental results have been obtained by Noir et al. 10 . One of their configurations cor- 
responds to a situation with a stable boundary layer. Their mean zonal flow measurements 
have been made at a fixed location z = 0.22 in a cylinder of aspect ratio a = 1.0686 for the 
parameters u = 1, e = 0.4, E = 10~ 4 . Their results are reported in figure 16(a) together 
with our computational results and Wang's prediction. 

One can see that numerical and experimental results are in good agreement but depart 
from Wang's curve. This is not surprising as for u = 1 and a = 1.0686, an inertial mode is 
resonantly excited whose structure can be seen in figure 16(c). In this figure, we can indeed 
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FIG. 16. Comparison between experimental, theoretical and numerical results in a cylinder of 
aspect ratio a = 1.0686 for e = 0.4, E = 10~ 4 , u = 1. (a) Mean zonal flow ug/e 2 at z = 0.22. Red 
points together with the red error bars are the experimental measurements of Noir et al. 10 using 
LDV technics, black solid line and blue dashed line are the numerical results and the theoretical 
prediction, respectively, (b) Contour of ug/e 2 in the (r, z) plane. The mean zonal flow is obtained 
by averaging the azimuthal velocity over three libration periods, (c) Instantaneous axial velocity 
contours (at time t = 2ir/u). 



observe a well-defined pattern in the bulk of the fluid which corresponds to the structure of 
the resonant inertial mode. The presence of this mode modifies the spatial structure of the 
mean zonal flow which now strongly depends on the axial coordinate, as illustrated in figure 

25 



16(b). 

In order to quantify this complex mean zonal flow we introduce the following global 
quantity: 

E c , 2 = ^eVfi 2 2 d\/ (26) 

which measures the mean kinetic energy per unit of mass associated with the zonal flow in 
the bulk. The volume V has been chosen such that it does not contain the viscous boundary 
layers. In particular, we have excluded from V the sidewall region where fi 2 strongly varies. 
Typically, for a cylinder of aspect ratio a = 2, we have taken boxes of size ranging from 
r x z = [0.5 x 0.6] to [0.8 x 0.9] to obtain a mean value of the mean kinetic energy and 
estimate error bars. The variation of the mean kinetic energy associated with the zonal flow 
as a function of the libration frequency u is reported in figure 17, where uncertainties are 
estimated from the variation of the mean value from one box to another. 



E, 



c,2 



1,2)(2,4) (2,6)(1,4) (1,6) 




FIG. 17. Dimensionless average kinetic energy E c ^ [defined by (26)] rescaled by e 2 as a function of 
the libration frequency u for E = 5.10 -5 , e = 0.1 and a cylinder of aspect ratio a = 2. The vertical 
dashed lines and the symbols (n, k) indicate frequency and wavenumbers of inertial modes (see 
table I). The dashed-dotted line corresponds to the predictions obtained from Wang's theory 11 . 



Figure 17 shows that E c ^ is peaked for particular values of u. We have been able to 
check that these particular frequencies are resonant frequencies for a cylinder of aspect ratio 
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n = 3 


k = 
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0.7586 


0.4370 


0.3052 


k = 
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1.2681 


0.8174 


0.5901 


k = 


3 


1.5518 


1.1152 


0.8406 


k = 


4 


1.7075 


1.3343 


1.0509 


k = 


5 


1.7975 


1.4916 


1.2222 


k = 


6 


1.8527 


1.6043 


1.3592 


k = 


7 


1.8886 


1.6860 


1.4680 


k = 


8 


1.9130 


1.7463 


1.5544 



TABLE I. Inviscid frequency of axisymmetric inertial modes as a function of their radial wavenum- 
ber n and axial wavenumber k for a cylinder of aspect ratio a = 2 (see Greenspan 1 ) . 

a = 2 as used in the simulation (see table I). The inviscid axisymmetric inertial modes are 
defined by two indexes n and k corresponding to radial and axial wavenumbers, respectively. 
Only the frequencies associated with an even number of half axial wavelength are observed. 
This is in agreement with the constraint imposed by the forcing. Indeed, libration induces 
an azimuthal flow which is even with respect to the mid-plane (z — 0). Modes with an odd 
number of half-axial wavelength have an odd azimuthal velocity component and therefore 
they cannot be forced directly. We can also notice that the peaks are the largest for the 
frequencies associated with the smallest wavenumbers and that frequencies corresponding to 
n > 3 are not visible. This is a viscous effect which tends to damp the largest wavenumber 
modes. Note finally that the general trend of E c ^ with respect to u still follows the theory 
by Wang 11 . 

Resonances can also be seen when the aspect ratio of the cylinder is varied for a fixed 
frequency. Figure 18 shows the variation of the zonal flow kinetic energy as a function of 
the aspect ratio for two different frequencies: u = 0.1 and oo = 1.27. For u = 0.1, E c2 
remains constant and equal to the prediction by Wang 11 as no inertial mode is excited for 
the parameters which have been considered. For oj = 1.27, a clear peak is observed at 
a = 2, in agreement with the excitation of the inertial mode of wavenumbers (n,k) = (1, 2) 
(see table I). The excitation of this mode is also responsible for the deviation from Wang's 
theory at a = 4. It is interesting to note that resonance could also decrease the zonal flow 
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kinetic energy as observed for a = 1.5 and a = 2.2 where modes (3, 4) and (2, 2) are excited 
respectively. 
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FIG. 18. Kinetic energy of the mean zonal flow for uj = 0.1 (red squares) and to = 1.27 (blue circles) 
for various values of the aspect ratio of the cylinder a = R/H and E = 5.10 -5 , e = 0.1. The red 
continuous line and the blue dashed line are respectively the theoretical value of the kinetic energy 
in the absence of inertial modes for to = 0.1 (E c ^ — 1.88 x 10~ 3 ) and w = 1.27 {E c $ ~ 2.23 x 10~ 3 ). 
The frequency to = 1.27 corresponds to the frequency of the mode (1,2) for a cylinder of aspect 
ratio a = 2 (see table I). 



We have seen that, when an inertial mode is excited, the zonal flow is no longer homoge- 
neous in the bulk [see figure 16(b)]. The radial structure of the angular velocity f2 2 of the 
zonal flow at a fixed axial location is illustrated in figure 19 when the mode (n, k) = (2, 1) 
is excited. In this figure, the effects of the variations of the Ekman number (figure 19(a)) 
and of the libration amplitude (figure 19(b)) are analysed. A clear dependence with respect 
to the Ekman number can be noted near the axis. This dependence is not associated with 
a viscous detuning of the inertial mode frequency as the viscous correction is negligible for 
the Ekman number we have considered 29 . Instead, we think that it could be due to strong 
nonlinear interactions occurring in the region close to the axis where the viscous shear layers 
emitted from the corners encounter. 

As the strength of the internal shear layers depends on the Ekman number, it is not 
surprising to observe an increase of the zonal flow when the Ekman number decreases. 
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FIG. 19. Angular velocity Q2 = u e/( r £ 2 ) of the zonal flow in a cylinder of aspect ratio a = 2 for 
lo = 1.27, (a) for e = 0.05 and various values of the Ekman number: E = 10 -4 (black continuous 
line), E = 8.10~ 5 (green dashed-dotted line), E = 6.10~ 5 (red dashed line), E = 4.10" 5 (blue 
dotted line); (b) for E = 8.10" 5 and various values of the amplitude of libration: e = 0.1 (black 
thick continuous line), e = 0.08 (green dashed-dotted line), e = 0.02 (blue dashed line), e = 0.01 
(red continuous line). 



Unfortunately, our range of accessible Ekman numbers did not allow us to provide a scaling 
for the zonal flow. Note that when the zonal flow is generated by tidal forcing, a scaling 
in E~ 3 / 10 was observed experimentally in a sphere 6 , and a scaling in E~ 3 ^ 2 was suggested 
numerically in a spherical shell 30 . 

In figure 19(b), we observe that the zonal flow remains proportional to e 2 even in the 
presence of inertial modes. The very small departure to the scaling in e 2 close to the axis 
where the amplitude is the largest is probably due to nonlinear effects. 
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VI. CONCLUSION 



In this paper, we have first analysed the occurrence of Taylor- Gortler vortices on the 
sidewall of a librating cylinder. We have shown that the vortices are generated by a cen- 
trifugal destabilization of the periodic boundary layer flow forced by the libration. We have 
derived a viscous stability criterion for the libration amplitude by considering the most un- 
stable azimuthal velocity profile obtained at leading order during a period near the sidewall. 
This criterion for the threshold amplitude provides a scaling in E 1 / 2 and a dependency 
with respect to u in qualitative good agreement with the numerical and experimental data. 
Quantitative differences have been observed and some explanations have been proposed. As 
the instability criterion only depends on the local structure of the boundary layer, which 
just has to be vertical and axisymmetric, it can be applied to other configurations. In par- 
ticular, we could apply the criterion to a librating sphere in the boundary layer close to 
the equator. The nonlinear regime of the instability has also been studied numerically. We 
have shown that the Taylor- Gortler vortices can force coherent structures in the bulk of the 
cylinder whose characteristics exhibit a preferred orientation around 45° independent of the 
libration frequency. 

Below the instability threshold, we have considered the mean zonal flow generated by 
libration. We have shown that the theory by Wang 11 captures the main characteristics 
of the zonal flow in the bulk for small frequencies and for frequencies larger than 2, i.e. 
when inertial modes are not directly excited. However, an unexplained discrepancy has 
been systematically observed close to the sidewall. The zonal flow in the presence of excited 
inertial modes has also been studied numerically. We have shown that the peaks of the 
zonal flow kinetic energy as the libration frequency or the cylinder aspect ratio is varied 
correspond to the resonant excitation of large-scale inertial modes. The dependence of the 
zonal flow with respect to the Ekman number and the libration amplitude has also been 
documented. 

All the numerical simulations have been carried out with an axisymmetric code. De- 
spite this constraint, we have seen that, in a librating cylinder, a mechanism of generation 
of inertial waves, different from direct forcing and period doubling, could be active when 
the sidewall boundary layer becomes turbulent. It would now be interesting to determine 
whether this process is active in 3D and how it depends on the geometry. In particular, a 
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precise description of the dynamics of librating ellipsoidal containers 31 ' 32 and the influence 
of the topography on the resulting flow would be of great interest. 
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APPENDIX 

The expression for the mean angular velocity, f^, appearing in formula (25) are given 
below for to < 2 and ou > 2 respectively (Note that these expressions were first derived by 
Wang 11 using a different notation). 

For u < 2: 

(u + 1) ^4 + 2 u + 2 (u - 1) V4 - 2 u - 2 1 

2 ~ 8 v / 4 + 2w(u; 2 + 4w + 5) ~ 8 ^4 - 2 w (w 2 - 4u + 5) ~ ^/I^Ju^/I+Ju 

V4-2o; V4 + 2w + w(v / 4 + 2w - ^4 - 2 lo) - 4 

+ V4-2w v /4 + 2w(2 + v / 4 + 2a; v / 4^ r 2"7) ( V4 + 2 w + v / 4 z: 27J) ' 
and for ou > 2 

+ 1) V2u + 4 + 2 (o;-3) V2a;-4 + 2 

2 ~ 8 v / 2w + 4(w 2 + 4w + 5) 8 V2 w - 4 (w 2 - 4u + 5) 

2u- ^2^ + 4^2^-4 2 w (w 2 - 6) - V2 w + 4 (w 2 + 3 w - 12) 

2 w v / 2aT+4 v / 2aT =r 4 4w(w 2 -3) V2W + 4^2^-4 

v / 2w-4(w 2 -3a;- 12) + 3u ^/2 u + 4 ^2 u - 4 
4w(w 2 - 3)^2^ + 4^277^4 ' 
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